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DEVELOPMENT OF A GLOBAL 3-D MAGNETOHYDRODYNAMIC 
COMPUTATIONAL MODEL FOR SOLAR WIND/COMETARY 
AND PLANETARY STUDIES 


1. INTRODUCTION 

This is the first Quarterly Progress Report under Contract 
NASW-4011. The work under this contract involves the development of a 
global, 3-D magnetohydrodynamic computational model to quantitatively 
describe the detailed continuum field and plasma interaction process of 
the solar wind with cometary and planetary magneto/ionopause shapes. 
The specific objectives of this first phase of the proposed three-phase 
investigation is to extend the present highly-successful solar 
wind/terrestrial planet interaction (Level 1) model (which is based on 
an axisymmetric gas dynamic plus frozen field approximation to the full 
MHD equations) to a full 3-D gas dynamic (Level 2) approximation, and 
also to develop and implement a mass loading capability in the Level 1 
interaction model . 

During this reporting period the technical tasks worked on included 
development and implementation of several enhancements to the present 
Level 1 model, preliminary development of the 3-D gas dynamic flow field 
solvers for the Level 2 model, and development of the 3-D frozen 
magnetic field solver for the Level 2 model. The work accomplished on 
each of these tasks is described in detail in the following sections. 


2. ENHANCEMENTS OF LEVEL 1 MODEL 

A number of important enhancements were developed and implemented 
in the Level 1 solar wind/terrestrial planet interaction model. First, 
the basic version of the Level 1 code that was developed on the 
NASA/Ames Research Center CDC 7600 computer facility was appropriately 
modified, brought up, and, through a series of case studies, verified on 
the RMA PRIME superminicomputer. A separate double-precision version of 
that code was then also developed and archived. This was done to enable 
exact comparisons to be made of results obtained from the Level 1 model 
when employing the RMA superminicomputer against corresponding model 
results obtained when using certain mainframe supercomputer facilities, 
such as the CRAY XM-P. For comparative calculations, a double-precision 
version is necessary due to the difference in significant digits carried 
for single precision arithmetic on CDC and CRAY supercomputers (14-16 
digits), and those on IBM and IBM derived operating systems (7-8 digits) 
such as the PRIME superminicomputer. Both single and double-precision 
versions of the Level 1 model are now fully active and, in fact, copies 
have been requested and distributed to several interested space 
scientists. Finally, an extended graphics capability was developed for 
the Level 1 model. This was done to provide a complete graphical output 
on the RMA Printronix printer of all the detailed plasma and field 
parameters calculated by the model throughout the interactive 


magnetosheath region. This task required an entire rewrite of the 
originally developed graphics. The output package now includes: maps 
of the plasma streamlines, velocity, temperature, and density contours; 
magnetic field contours and field line patterns; spacecraft trajectory 
diagrams (ecliptic and polar plane views); time histories of plasma 
density, temperature, velocity magnitude, three velocity components; 
magnetic field magnitude, and three magnetic field components. This 
feature of having a complete graphics output package for use on a 
high-density dot matrix printer, has never been available before with 
the model. The output is essentially report quality, and the capability 
provides in effect the means for instantaneous viewing of the model 
results at almost no cost. This is ideal for comparative studies with 
other theories and, particularly, observational results. Figures 1-3 
provide some examples of an abbreviated subset of the complete graphical 
output. Finally, the graphics package developed was intentionally based 
on CALCOMP routines which is a widely available industry standard. This 
enables the entire package to be easily transferable to other 
facilities. 


3. PRELIMINARY DEVELOPMENT OF LEVEL 2 FLOW FIELD MODEL 

For reasons that were described in detail in the proposal, the 
general procedure for the development of the model under this contract 
is such that in its final form the model will embody a hierarchy of 
three different levels of accuracy and associated computational 
algorithms. These levels are: 


Level 


Model 


1 Ax i symmetric gas dynamic + frozen B field 

2 3-D gas dynamic + 3-D frozen B field 

3 3-D MUD 


A major objective of the first phase of the present study is to develop 
the Level 2 3-D gas dynamic plus 3-D frozen B field model. For reasons 
also described in detail in the proposal, all 3 levels of the model will 
contain two separate flow field solvers coupled together to determine 
the entire steady state flow field. This feature of our model, which is 
unique among all currently existing global interaction models, allows 
very high resolution of the entire flow field, particularly in the tail 
region. As illustrated in Figure 4, where a typical flow field grid 
density of the combined Level 1 axisymmetric flow solvers is shown, the 
region from the subsolar poirt to the terminator plane is solved via a 
time marching solver (NOSE code) to the steady state; whereas the 
region downstream of the terminator is solved via a spatially marching 
solver (TAIL code) which advances the solution downstream as far as 
required. The development of the corresponding 3-D versions of the NOSE 
and TAIL codes (N0SE3D, TAIL3D) is described below. 



3.1 Level 2 Nose Region Flow Solver: N0SE3D Code 


After an extensive review, the selection of the basic method to be 
employed in the N0SE3D code has now been made. The methodology will be 
based upon the AIR3D code developed at NASA/Ames Research Center 
(Ref. 1). The solution algorithm employed is the Beam and Warming 
implicit approximate factorization algorithm used in delta form 
described in Ref. 2. The original basic algorithm is first-order 
accurate in time, noniterative, second-order accurate in spatial 
derivatives with central difference operators. This methodology is 
capable of solving both the inviscid gas dynamic Euler equations as well 
as, if required, a certain form of the viscous gas dynamic equations 
known as the thin-layer Navier-Stokes equations. These equations 
include certain additional terms over and above the Euler equation level 
which represent viscous effects in the body-normal direction. This 
capability is important to our study for future planned applications in 
order to study certain viscous effects near magneto/ionopause 
boundaries. The mesh employed in this method is one which incorporates 
a fitted shock outer boundary which exactly satisfies the unsteady, 3-D 
Rankine-Hugoniot relations and moves with the shock during the 
convergence process. The inner boundary is implemented as a fitted 
stationary surface. A mesh clustering capability for high density 
clustering in the body-normal direction has been implemented 
specifically to enable accuracy to be maintained for viscous 
calculations. The N0SE3D code is currently being modified for 
implementation on the RMA PRIME computer. Additionally, comparative 
computations are underway using the Level 1 NOSE code on certain simple 
magneto/ionopause shapes. These are being archived for use in future 
benchmark comparisons with the Level 2 N0SE3D code. 


3.2 Level 2 Tail Region Flow Solver: TAIL3D Code 

The gas dynamic solver to be employed for the Level 2 solution of 
the 3-D tail region of the flow field will be based on the 3-D gas 
dynamic marching code reported in Ref. 3. The solution algorithm is a 
shock capturing, fully conservative form that employs second-order 
noncentered spatial derivatives to solve the steady, inviscid, gas 
dynamic Euler equations. As with the N0SE3D code, the mesh to be used 
with the T4IL3D code will employ a fitted outer boundary for the bow 
shock, at which location the steady 3-D Rankine-Hugoniot relations are 
exactly satisfied. The inner boundary of the mesh will also be fitted 
but to an impenetrable surface. Currently, the necessary geometry 
routines required to apply the TAIL3D solver to nonaxisymmetric 
magneto/ionopause surfaci_s typical of those in the solar system (Earth, 
Jupiter, Saturn) are being developed. 


3.3 Development of Level 2 3-D Frozen B Field: MAG3D Code 

A method has also now been chosen to determine the frozen 3 field 
for a general 3-D gas dynamic flow field. The implementation involves 
the 3-D generalization of the kinematic procedure used in the present 


Level 1 model to determine the perpendicular component of the 
magnetosheath £) field via the Alksne and Webster (Ref. 4) decomposition 
method. That simpler decomposition method is valid only for 
axisymmetric flow fields and no longer applies to 3-D flows. The 
generalized 3-D kinematic g field procedure has now been derived. The 
coding of a preliminary version of the method is almost complete. 


4. TASKS PLANNED FOR NEXT REPORTING PERIOD 

During the next reporting period, work will be concentrated on 
completion of the 3-D frozen field (MAG3D) code. Comparative 
computations of that generalized procedure will be made with the 
decomposition method in the present Level 1 model. Modifications will 
be completed on the 3-D nose region flow solver (N0SE3D) to enable the 
code to be brought up on the RMA PRIME computer. The code will then be 
activated and initial verification calculations carried out as far as 
possible. Separate benchmark computations of the N0SE3D code will be 
also carried out on the NASA/ARC CRAY XM-P for future comparative 
studies. Finally, the special geometry routines needed to describe 
nonaxisymmetric magneto/ ionopause shapes will be developed for the 
TAIL3D code. If time permits, the TAIL3D code will be brought up on the 
RMA PRIME computer. 
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Figure 2. Sample magnetic field output from Level 1 model graphical output package: 
magnetic field lines and contours for perpendicular field component for 
fi, = 3.0, y = 5/3, H/R 0 = 0.03 ionopause obstacle flow 
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Figure 3. Sample time history output from Level 1 model graphical 
output package: temperature vs. time for a trajectory 

passing through the magnetosheath region and ionopause 
of the flow field shown in Figure 1 



Figure 4 . Illustration of typical flow-field grid density for 
gas dynamic solution; H ro = 3.0, y = 5/3, H/R 0 = 0.03 



